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Simple  Frequency  Estimation 
via  Exponential  Samples 

Steven  Kay,  Fellow,  IEEE 


Abstract — A  method  for  determining  the  frequency  of  a  real 
sinusoid  is  proposed.  Based  on  exponential  sampling  of  the 
waveform,  the  approach  requires  virtually  no  computation.  It  can 
be  easily  implemented  in  digital  hardware. 


I.  Description  of  Method 


ASSUME  we  observe  the  real  continuous-time  sinusoid 
s(t)  =  A  sin  (2nFot),  where  A  >  0  and  that  we  wish  to 
determine  the  frequency  Fo,  which  satisfies  0  <  Fq  <  Fmax- 
The  sinusoidal  phase  is  assumed  to  be  zero.  Let  Fmax  =  2M 
for  some  integer  M.  To  measure  the  frequency  to  an  accuracy 
of  N,  we  use 


where 


N 

Ft.  =  iW  X>n2-n 

n= 1 

bn  =  0  if  s(fn)  >  0 

=  1  if  s(i„)  <  0 


(1) 


and  tn  =  2"  M  l.  If  s(tn)  =  0  for  some  n  =  no,  then 


/n°- 1 

-Fo  =  Fmax  (  bn2~n  +  2' 


■n0 


„  n=l 


(2) 


which  effectively  terminates  the  expansion. 


II.  Rationale 

Define  the  normalized  frequency  as  /0  =  (Fo/Fmax)  so 
that  0  <  f0  <  1.  We  use  a  binary  expansion  for  /0  as 

OO 

f0=Y/bn2~n  (3) 

n=l 


angle  of  the  sinusoid  or  9{t)  =  27rFof.  For  sample  times 
tn  =  2n_M_1,  the  angle  samples  become 

0n  -  27rFofn 

=  27r/0FMAX2”-M-1 
=  27r/o2n_1  (4) 


for  n  =  1,2,---  .  Consider  n  =  1  so  that  0 j  =  2nfo ■  If 
6i  =  0,  then  from  (3),  0  <  /o  <  1/2.  The  upper  limit 
of  1/2  is  assumed  not  to  be  possible  due  to  assumption 
of  a  terminating  representation  since  1/2  is  represented  as 
0.1000  •  •  •  .  Hence,  &i  =  0  if  and  only  if  0  <  9i  <  7r.  In 
addition,  if  b\  =  1, 1/2  <  /o  <  1  so  that  6j  =  1  if  and  only 
if  it  <  9 1  <  2it.  Next,  consider  n  —  2  so  that  9 1  —  Aitfo. 
Using  (3) 


=  4,(|  +  |>2- 

=  2  tcbi  +  2tt  ^  bn2~n+1 


We  can  omit  the  27r6j  term  since  it  equates  to  0  or  27 r,  and 
thus,  it  does  not  affect  the  sine  function.  Let  a  =  b2/2  + 
££L36n 2~n+1  and  note  that  if  b2  =  0,0  <  g  <  1/2  and  if 
b2  =  1,1/2  <  a  <  1.  Hence 


0  <  92  <  it  if  and  only  if  b2  =  0 
tt  <02  <  2ir  if  and  only  if  b2  —  1. 


By  continuing  the  same  argument,  we  can  show  that 


0<#7i<7T4=>bn=0 
it  <  9n  <  2-k  bn  ~  l.  (5) 


Excepting  the  case  when  9n  =  0  or  7r,  this  translates  into 


where  bn  =  0  or  1. 

For  uniqueness,  it  is  assumed  that  the  dyadic  rationals, 
that  is,  numbers  of  the  form  k/2h  for  L  an  integer  and 
k  =  0, 1,  •  •  ■ ,  2l  —  1,  are  represented  by  terminating  expan¬ 
sions.  Hence,  the  number  3/4  is  represented  as  0.11000  •  -  - 
as  opposed  to  0.10111  •••  .  With  (3),  determination  of  /o 
is  equivalent  to  determination  of  {bn}.  Now,  consider  the 
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s(tn)  >  0  44-  bn  =  0 
s(tn)  <  0  O’  bn  =  1. 

If  s(f„0)  =  0  for  some  n  =  no,  then  9no  =  irk  for  k  an 
integer.  Thus,  from  (4) 

irk  —  27r/o2"0_1 

or  fo  =  k/2n°,  which  is  a  dyadic  rational.  It  then  follows  that 
9n  =  2irf02n-1 

=  2it^2n-x 
2n° 

=  irk2n~n°  (6) 
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which  is  a  multiple  of  27t  for  n  >  n0,  and  hence,  6n  =  0  for 
n  >  n o  based  on  (5).  If  no  is  the  sample  for  which  s{tn)  is 
first  equal  to  zero,  then  k  must  be  odd.  Otherwise,  k  =  21, 
and  from  (6) 

9n  =  7rZ2n-no+1 

so  that  s(f„0_i)  =  0,  contradicting  the  assumption  of  the  first 
zero.  In  summary,  if  no  is  the  first  sample  term  for  which 
s(f„)  =  0,  then  we  should  choose  bna  =  1  since  6no  =  bn 
for  k  odd,  and  bn  =  0  for  n  >  no  since  9n  will  be  a  multiple 
of  27T. 

As  an  example,  assume  Fmax  =  64  Hz  and  Fo  =  60  Hz 
so  that  M  =  6. 

For  N  =  4  bit  accuracy,  we  sample  at 

tn=  2n~M~1  n  =  1,2, 3,4 

_  i_  J_  _L  l  Q 
—  64  >  32  >  16’  8 

Then,  s(fn)  =  Asin(27r(60)fn)  so  that  s(fi)  <  0,s(t2)  < 
0,  sfa)  <  0,  s(f4)  =  0.  Applying  (2)  with  no  =  4,  we  have 

F0  =  64(2-1  +  2-2  +  2-3  +  2~4)  =  60  Hz. 

III.  Practical  Considerations 
A.  NoiselPhase  Errors 

Assume  that  the  sinusoidal  phase  is  not  zero  so  that  s(t)  = 
A  sin  (2tt F0t+4)).  We  model  $  as  a  random  variable  uniformly 
distributed  over  [—a,  a]  or  0  ~  U[—a,  a].  In  addition,  due  to 
observation  noise,  we  observe  the  samples 

x(tn)  =  Asin(27rFof„  +  <j>)  +  w[n] 

where  io[n]  is  modeled  as  white  Gaussian  noise  with  variance 
a1  and  is  independent  of  <j>.  The  mean  squared  error  (MSE) 
of  Fo  as  given  by  (1)  is  now  found.  We  need  not  consider  (2) 
since  x{tn)  /  0  with  probability  one.  The  mean  squared  error 
of  F0  =  FMax  £n=i  b„2~n  is 

MSE(F0)=F[(Fo-Fo)2] 

=  EtEwl4(Fo-F0)2} 

=  E<t>Ew{(F0-F0)2} 

where  E#  denotes  the  expected  value  with  respect  to  the  PDF 
of  <f> ,  and  Ew  denotes  the  expectation  with  respect  to  the  w[n] 
samples.  However 

EW[(F0  -  Fo)2]  =  var^Fb)  +  (EW(F0)  -  F0)2. 

Since  bn  will  be  1  when  x(tn)  <  0  and  0  otherwise,  we  have 
that  Pr[6n  =  1]  =  Pr[s(t„)  <  0]  =  Q((A/cr)sin(27rFofn  + 
0))  =  Pn{4>)  and  therefore,  conditioned  on  cj>,bn  is 

Bernoulli  with  probability  pn(<f>).Q(u)  is  defined  as  Q(u)  = 
/“  (l/x/2^)e-^1^t2  dt.  Hence 

N 

var ^(Fo)  =  F^AX  ^  varu,(6„)2-2n 

n=l 

N 

=^AxEp«w(1-p«w)^2,, 

n=l 


U.6,  N-8,  FO-59 


16 

SNR  In  dB 


Fig.  1.  MSE  for  a  frequency  of  Fo  =  59  Hz. 
and 

N 

EW(F0)  =  Fmax  y~^Pn(^)2~n. 

n=l 

Since,  <fi  ~  U[—a,a],  we  finally  have  that 

MSE(Fo)  =  Ya  fa  fMA xEp»  -  Pn(m-2n  # 

+  ^/0  (irMAX^PnW2-n-Fo^  d<j> 


where 


=  Q(jsi 


sin  (27T F0t„  +  <f>)  . 


As  an  example,  the  MSE  is  shown  in  Fig.  1  for  a  frequency 
of  Fo  =  59  Hz.  The  maximum  frequency  was  chosen  to  be 
Fmax  =  64  Hz  (M  =  6)  and  N  =  8  bits  of  resolution  were 
selected. 

B.  Observation  Interval 

Since  the  sample  times  are  exponential,  the  length  of 
time  over  which  we  must  observe  x(t)  can  become  large. 
The  observation  interval  for  N  bits  is  T  =  — 

2n~  1  /Fmax  A  given  relative  bias  error  due  to  bit  truncation 
or  |F(F0)  -  Fo|/F0  =  e  may  be  obtained  if  we  choose  N 
bits  where 

2-^Fmax  _ 

F0  _£' 

(This  ignores  the  bias  due  to  noise.)  Since  2T  =  2n/Fmax, 
we  have 


The  observation  time  is  determined  by  the  accuracy  desired 
and  the  minimum  frequency  to  be  estimated. 


KAY:  SIMPLE  FREQUENCY  ESTIMATION 

rv.  Discussion 

The  approach  described  is  easily  implemented  in  digital 
hardware  using  a  clock  running  at  a  rate  of  Fmax  cycles/s  and 
a  binary  counter  to  obtain  the  exponential  sample  times.  The 
samples  are  taken  when  the  counter  contains  1,  2,  4,  8,  •  ■  •  . 
To  ensure  a  sinusoidal  phase  of  zero,  the  counter  is  initialized 
when  x[t )  crosses  zero  in  an  upwards  direction.  From  (2), 
the  probability  of  x(tn)  =  0  and,  therefore,  of  terminating 
the  expansion  is  zero  in  the  presence  of  noise.  Thus,  if  f0  is 
a  dyadic  rational,  whose  binary  expansion  should  terminate, 
noise  effects  may  cause  the  expansion  to  continue.  To  avoid 
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this  it  may  be  better  to  truncate  the  expansion  at  n  =  n0  if 
|a;(fno)|  <  6,  where  8  is  a  small  number. 

Last,  the  estimation  method  may  be  used  to  measure  the 
frequency  of  other  periodic  waveforms  such  as  square  waves, 
for  example.  The  only  requirement  is  that  we  should  be  able 
to  determine  from  the  waveform  if  9n  takes  on  values  as  per 
(5).  In  addition,  the  speed  of  a  motor  or  rate  of  revolution  of 
any  constant  rotation  object  may  be  ascertained  as  described. 

Acknowledgment 

The  author  wishes  to  thank  S.  Talwalkar  for  providing  the 
computer  results. 


Methods  for  Chaotic  Signal  Estimation 


Steven  Kay 

University  of  Rhode  Island 
Kingston,  RI  02881  * 

Venkatesh  Nagesha 

Speech  Research  Department 
AT&T  Bell  Laboratories 
Murray  Hill,  NJ  07974 

Abstract 

A  dynamic  programming  algorithm  and  a  suboptimal  but  computation¬ 
ally  efficient  method  for  estimation  of  a  chaotic  signal  in  white  Gaussian 
noise  are  proposed.  The  nonlinear  map  is  assumed  known  so  that  only  the 
initial  condition  need  be  estimated.  Computer  simulations  confirm  that  both 
approaches  produces  efficient  estimates  at  high  signal-to-noise  ratios. 

1  Problem  Statement 

The  problem  is  to  estimate  a  chaotic  signal  s[n]  that  is  embedded  in  white  Gaussian 
noise  io[n]  of  variance  a2.  Hence,  the  data  model  is 

rcjn]  =  s[n]  +  w[n]  n  =  0, 1, . . . ,  N  —  1,  (1) 

where  s[n]  is  a  deterministic  signal  generated  according  to  the  nonlinear  and  non- 
invertible  map 

SM  =  f(s[n-  1])  .  (2) 

EDICS  #  2.7.5,  3.6.1 
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The  map  /  :  [0,  1]  — >  [0,  1]  is  assumed  known  and  to  be  unimodal.  As  a  result,  it 
has  two  preimages.  Typical  examples  of  such  maps  are  the  logistic  and  tent  maps 

[1] .  We  furthermore  assume  that  the  initial  condition  s[0]  is  unknown  so  that 
the  problem  of  estimation  of  s[n]  reduces  to  that  of  estimation  of  only  the  initial 
condition.  Once  this  estimate  is  available  the  entire  signal  may  be  estimated,  at 
least  in  theory,  using  (2).  In  practice,  because  of  the  extreme  sensitivity  to  errors, 

(2)  cannot  be  used  to  propagate  signals  in  the  forward  direction.  Instead,  we  rely  on 
the  parameterization  described  in  [2]  whereby  the  sequence  {s[0],  s[l], . . . ,  s[N—  1]} 
is  replaced  by  {po,pi, . . . , p/v-2,  s[N  —  1]}.  The  sequence  {pn}n=o  is  the  itinerary, 


0  if  0  <  s[n]  <  c 

Pn=  {  , 

1  if  c  <  s[n ]  <  1 

where  c  is  the  value  of  x  for  which  /(x)  is  maximum.  Hence,  the  itinerary  can  be 
used  to  determine  the  appropriate  preimages  of  a  chaotic  signal  when  propagated 
backward.  As  an  example,  for  the  logistic  map  or  /( x)  =  4x  (1  —  x),  the  function 


i  +  (2p- i)yr-x 

fP  ix)  = - 2 - 


and  consists  of  the  preimages  f0  1  (x)  = 


1— \/l—  X 


and  fx  J(x)  =  Thus, 


can  write 


s[n  —  1]  =  /^(sM)  . 


2  Halving  Method 

Before  discussing  the  dynamic  programming  (DP)  approach  to  the  maximum  like¬ 
lihood  estimate  (MLE)  evaluation,  we  describe  a  suboptimal  but  computationally 
efficient  approach.  In  the  process  we  prove  that  for  a  large  class  of  unimodal  maps 
the  initial  condition  can  be  determined  from  the  itinerary  of  the  map.  The  ar¬ 
gument  relies  on  certain  topological  conjugacy  relations.  Now  assume  that  /  is 


2 


0  <  x  <  | 

i<*<i 

so  that,  there  exists  an  invertible  transformation  h  :  [0,  1]  — *  [0,  1]  such  that 

f  =  hoToh~K  (3) 

Then,  we  may  write  s[n]  as 

s[n]  =  /n(s[  0]) 

=  hoTn  o  A_1(s[0])  , 

since  fn  =  h  o  Tn  o  h -1  and  /"  is  the  n  fold  composition  of  /.  Let  s'jO]  =  h~x  (s[0]) 
so  that, 

Tn(s'[ 0])  =  h-\s[n})  . 

Now  use  the  substitution  property  or  Tn  =  T  o  5'n~1  [3],  where  S  is  the  Bernoulli 
shift 

S(x)  —  2smod  1  0  <  x  <  1  , 

to  yield 

T  o  5'n_1(s,[0])  =  /i-1(s[n])  . 

If  s'[0]  is  represented  in  the  binary  format  s'[0]  =  ^*2  k  =  0.bib2  . . .,  (where 

we  use  a  terminating  expansion  for  the  dyadic  rationals),  we  have  that 

Sn-\S,[Q])=  O.bnbn+l...  ■ 

Hence,  it  follows  that 

T(0.6n6n+1 ...)  = /i_1(s[n])  . 

The  effect  of  the  tent  map  is  to  shift  left  by  one  bit  if  the  argument  is0<x<l/2 
and  shift  left  and  complement  (due  to  1  —  x  factor)  if  1/2  <  x  <  1.  As  a  result, 

^-1(s[n])  =  0.6n+i&n+2  . . .  if  bn  =  0, 

=  0.6n+l^n+2  •  •  •  if  bn  —  1) 


conjugate  to  the  tent  map  T  or 


T(x)  = 


2x 

2(1  -  x) 


3 


where  the  overbar  denotes  the  complement  bit.  This  says  that 

f 

0 

if  0  <  /i-1(s[n])  <  | 

and 

bn  =  0 

&n+l  —  " 

1 

if  |  <  /i-1(s[n])  <  1 

and 

o 

II 

-o 

1 

if  0  <  /i-1(s[n])  <  \ 

and 

bn  =  l 

0 

if  |  <  /i_1(s[n])  <  1 

and 

bn  =  1 

If  we  now  let  pn  be  the  itinerary  of  h  1(s[n])  or  pn  =  0  if  0  <  h  1(s[n])  <  1/2 
and  pn  =  1  if  1/2  <  /i-1(.s[n])  <  1,  we  have  finally  that, 

bn+1  =  ®  Pn  i 

where  ©  denotes  the  exclusive  OR  operation.  The  recursion  begins  with  bx  ~  p0. 
To  summarize,  we  can  determine  the  initial  condition  of  a  map  that  satisfies  the 
conjugacy  relation  of  (3)  as 

s[0]  =  />(f>n2— 

\n=l 

where  bn  =  6n_x  ©  pn_i  ,  ( b\  =  Po)  and  pn  is  the  itinerary  of  /i-1(s[n]).  It  is 
interesting  to  note  that  in  practice  we  will  have  {p0,pi, . . .  ,pw-i}  based  on  the 
itinerary  of  the  data  set  {s[0],s[l], . . .  ,s[N  —  1]},  assuming  no  noise.  Hence,  we 
will  have 

so  that  the  estimate  of  /i-1(s[0])  will  be  in  error  by  at  most  l/2N .  We  will  call 
the  estimation  of  s[0]  by  (4)  the  halving  method  since  as  each  itinerary  value  is 
obtained,  the  interval  in  which  s'[0]  must  reside  is  halved. 

As  an  example,  for  the  chaotic  logistic  map  /(x)  =  4x(l  —  x )  it  is  known 
that  h(x)  =  sin2(:|x)  and  h~1(x)  =  ^arccos(l  —  2x).  Since  h~l  maps  [0, 1/2)  onto 
[0, 1/2)  and  [1/2, 1)  onto  [1/2, 1)  the  itineraries  of  /i._1(s[n])  and  s[n]  are  the  same. 
Thus,  the  halving  method  estimates  s[0]  as 

s[0]  =  sin2  ,  (5) 

A  A 

where  bn  =  0  pn-\  and  pn  is  the  itinerary  of  x[n].  The  performance  for 

this  example  will  be  discussed  in  Section  4.  Note  that,  if  we  consider  s[n]  as  the 
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initial  iterate,  then  it  depends  only  on  {pn,pn+i, . . Furthermore,  s[n]  may  be 
determined  directly  from  this  itinerary  in  a  similar  fashion  as  per  (4).  Or  since 
the  itinerary  {pn,Pn+i,  •  •  •}  can  be  obtained  for  any  s[n]  there  is  a  one-to-one 
correspondence  between  s[0]  and  {po,Pi,  ■  •  •  ,Pn-2,  ■s[Ar  —  1]}  as  already  noted  in 
[2]- 

3  Dynamical  Programming 

The  MLE  for  the  initial  condition  is  obtained  as  the  value  of  s[0]  that  minimizes 
/  \2 

J  =  2_,  (  x[n]  —  s[n]  ]  where  s[n]  =  /n(s[0]).  A  straightforward  minimization 

n= 0  '  ' 

requires  one  to  compute  /n(.s[0]),  which  leads  to  computational  errors.  Rather, 
we  employ  a  DP  approach,  which  does  not  require  a  forward  propagation.  The 
method  to  be  described  applies  to  any  unimodal  map.  In  particular,  for  the  tent 
map,  which  is  piece-wise  linear,  an  analytical  solution  may  be  found  as  in  [2].  Let 


and  also 


Jk  (#])  =  J2  (xN  -  5N)2 


s\n]  =  /  1  r  (sfArl) 

Pn  iPn— 1  1  '  *■  *' 


be  the  inverse  function  composition  for  n  <  k  —  1 .  This  is  defined  as 

s[k-l]  =  /"!,(#]) 
s[k-2]  =  f-?_,{s[k- 1]) 

=  fpLifpLMW) 


Now,  let 


Jfc(s[A:])  =  min  Jk(s[k]) 

{PO,Pl,--,Pk-l} 


so  that  the  desired  minimization  is  effected  when  k  =  N  —  1  and  h(s[N  —  1])  is 
minimized  over  s[N  —  1].  Thus, 

h{s[k])  =  min  Y  (XN  “  5W)2 

iPo,Pi,-,Pk-i)  n=0 
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k-1 

=  min  min  (z[nj  —  .s[ra])2  +  ( x\k ]  —  s\k])2 

Pk-l  {po.Pli—,Pfc- 2}  n=0 

=  min  min  [jfc_i(s[fc  —  1])  +  (x[fc]  —  s[£])2j 

Pk-l  {po,Pl,-;Pk-2}  *■  J 

Since  s[k]  does  not  depend  on  {p0,Pi, . . .  ,pk-i}  but  only  on  {pk,...,PN-2,s[N  — 
1]},  as  per  (6)  we  have, 

Ik(s[k])  =  min Ik-i  (.££.,(«[*]))  +  (*[*]  -  «[*])*  ,  (7) 

as  our  DP  algorithm.  The  recursion  is  computed  as  a  function  of  s[k]  for  k  — 
1,2,...,  TV  —  1  with  the  initialization 

'  lo(s[0])  =  (m  -  s[0})2  . 

After  obtaining  In-i{s[N  —  1]),  we  find  the  s[N  —  1]  that  minimizes  it  and  denote 
the  minimizing  value  as  s[N  —  1].  We  then  backtrack  to  find  {po,Pi>  •  •  • ,  2} *  In 

doing  so,  the  signal  sequence  is  estimated  as 

s[n-l]  =  /r1_i(s[n]) 

for  n  —  N  —  1,  N  —  2, . . . ,  1.  Note  also  that  this  approach  avoids  the  exponential 
increase  in  computational  errors  since  the  propagation  is  along  the  stable  manifold. 

4  Computer  Simulations  Results 

It  has  been  shown  [4]  that  the  MLE  of  the  initial  condition  of  a  one-dimensional 
chaotic  signal  is  an  inconsistent  estimator.  Hence,  the  usual  asymptotic  distribu¬ 
tion  as  N  — >  00  does  not  hold.  However,  as  the  SNR  becomes  large,  the  asymptotic 
distribution  is  valid.  In  particular,  the  asymptotic  probability  density  function 
(PDF)  of  s[0]  is 

l[0]~V(S[0],7-1MO])),  (8) 

where  7(s[0])  is  the  Fisher  information  (and  hence,  J-1(s[0])  is  the  Cramer- Rao 
lower  bound  (CRLB)  for  an  unbiased  estimator)  and  is  found  to  be 

^(-[OD  =  >  (9) 

Zjm=0  Pm 
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where 


* 

1  if  m  =  0 

<  771  —  1 

n  i/'  (SW)]2  if  m  >  ° 

.  /=0 


and  f  denotes  the  derivative  of  /  [5]. 

We  compare  the  performance  of  an  MLE  using  a  fine  grid  search  (by  minimizing 
iV-l 

J=  E  (x[n]  —  s[n])2  directly)  and  that  of  the  two  algorithms  previously  described 

n=0 

to  the  theoretical  performance  of  (8).  In  particular,  the  bias  and  variance  are 
determined  using  Monte  Carlo  simulations.  Note  that  the  CRLB  for  a  logistic 
map  can  be  shown  to  be 


r\*[  o])  = 


<725[1] 


We  used  a  data  record  length  of  N  —  10  and  an  initial  condition  of  s[0]  =  0.61. 
Some  implementation  details  for  the  three  methods  are  now  discussed.  For  the 
grid  search  MLE  we  first  use  a  coarse  search  by  dividing  up  the  [0, 1]  interval 
into  1000  points.  Then  we  search  over  the  interval  [sc[0]  —  0.001,  sc[0]  -f  0.001] 
using  1000  points  where  sc[0]  is  the  coarsely  obtained  point  of  the  minimum.  The 
implementation  of  the  halving  method  departs  slightly  from  (5).  This  is  because 
the  use  of  (5)  may  produce  a  maximum  error  of  1/2N  «  0.001,  even  in  the  absence 
of  noise.  The  CRLB,  because  of  its  exponential  decrease  with  N,  can  be  much  lower 
than  this  error.  To  improve  the  performance  we  first  obtain  the  coarse  estimate 


340]  =  sin2  6„2-“  +  2-("+'))) 


where  the  addition  of  the  term  2  (JV+1)  has  the  effect  of  choosing  the  midpoint 
/  n  N  j  \ 

of  the  interval  Then,  a  fine  search  (100  points) 

\n=l  n=l  ^  ) 

over  the  interval  (sc[0]  —  2-iV,sc[0]  +  2~N )  is  made  by  minimizing  J,  as  in  the 
MLE  implementation.  Note  that  the  length  is  twice  that  of  the  coarsely  obtained 
interval.  The  fine  search  is  repeated  twice  over  successively  smaller  intervals  of  100 
points.  The  DP  implementation  of  the  MLE  uses  500  points  for  the  state  s[k].  For 
the  search  over  s[N  —  1]  we  use  an  initial  quantization  of  500  points  to  locate  a 
coarse  estimate  of  s[N  —  1].  We  then  choose  a  finer  grid  for  s[iV  —  1]  centered  about 
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the  previous  estimate  as  well  as  the  same  estimates  {po,pi, . .  • , PN-2}  obtained 
previously. 

The  results  are  shown  in  Figure  1.  Above  the  threshold  of  about  40dB  all 
three  methods  produce  unbiased  estimates  that  achieve  the  CRLB.  (The  bias, 
although  not  shown,  was  negligible  above  the  threshold).  This  is  in  accordance 
with  the  theoretical  results  of  [4].  Of  the  three  methods  the  halving  approach  is 
computationally  the  simplest.  If  enough  data  points  are  available,  the  fine  search 
employed  for  the  halving  method  may  be  eliminated.  This  is  because  the  error  will 
be  at  most  1/2N  and  so  the  estimate,  which  although  not  attaining  the  CRLB, 
will  be  accurate  enough  for  most  practical  purposes. 
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